#ifndef LAB2_CRT_HELPER 
#define LAB2_CRT_HELPER 





template<typename T>
T gcd(const T& lhs, const T& rhs) 
{
    T n(lhs),m(rhs);
    if(n == T(0))
        if(m != T(0)) return m;
        else throw "There is no gdc between two zeroes";

    if(m == T(0)) return n;
    size_t d(0);
    while (((m & T(1)) == T(0))
        && ((n & T(1)) == T(0))) 
    { m >>= 1; n >>= 1; ++d; }

    while ((m & T(1)) == T(0)) m >>= 1;
    while ((n & T(1)) == T(0)) n >>= 1;

    while (m != n)
    {
        if (m < n) 
        {
            n -= m;
            do  n >>= 1;  while ((n & T(1)) == T(0));
        } 
        else
        {
            m -= n;
            do  m >>= 1;  while ((m & T(1)) == T(0));
        } 
    }
    return (m <<= d);
}


template<typename T>
T inv(T x, T mod)
{
    if ((gcd(x, mod) != T(1)) || (mod==T(1))) return T(0);
    T q, r, d, c;
    T  a, a1 /*, b, b1 */; 
    d = x % mod, c = mod;

    a = /* b1 = */ T(1), a1 = /* b = */	T(0);
    while(true)
    {
        q = c / d, r = c % d;
            
        if(r == T(0)) return a;

        c = d, d = r;

        T t(a1);
		a1 = a,
		a =  (t - q * a) % mod;
		
		if(a < 0) a += mod;
        //t = b1, b1 = b, b = t - q * b;
    }
}

#include <vector>
#include <numeric>

#include "..\Lab1_Large_numbers\__uintX2.h"

template<typename T>
T solve_congruences_system_with_CRT(std::vector<T> rems, std::vector<T> mods)
{
	const size_t congruences_number = mods.size();

	if(congruences_number != rems.size()) 
		throw "The number of remainders should be the same as the number of modules";

	__uintX2<T> M(1);
	for(size_t i = 0; i < congruences_number; ++i)
	{	
		M *= mods[i];
	}
	/*__uintX2<T> M =  std::accumulate(mods.begin(), 
			               mods.end(),
						   (__uintX2<T>)1, 
						   std::multiplies<T>());
						   */

	__uintX2<T> S(0);
	/*
	printf("\nM = %u\n",M); 
    printf("\n m\t M'\t M'^(-1)\t S\t rem\t"); 
    */
	for(size_t i = 0; i < congruences_number; ++i)
	{	
		//       n               -1
		// S = sigma  a * M' * (M'  mod m )
		//      i=1    i   i     i       i
    
		__uintX2<T> _M = M / mods[i]; 

		
		T _M_inv = 
			_M.inv(mods[i])._l
			//inv(_M , mods[i])
			;

		S +=  __uintX2<T>::mult_m( (_M * rems[i]), _M_inv, M);   
		
		
		//printf("\n %u\t", mods[i]); 
		//printf(" %u\t", _M); 
		//printf("M inv is %u\t", _M_inv); 
		//printf(" %u\t", rems[i] * _M * _M_inv); 
		//printf("%u\t",  (_M * _M_inv) % mods[i]); 
		
	}

	return (S % M)._l;
}
#endif /* LAB2_CRT_HELPER */ 
